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Quantum algorithms are well-suited to calculate estimates of the energy spectra for spin lattice 
systems. These algorithms are based on the efficient calculation of discrete Fourier components of 
the density of states. The efficiency of these algorithms in calculating the free energy per spin of 
general spin lattices to bounded error is examined. We find that the number of Fourier components 
required to bound the error in the free energy due to the broadening of the density of states scales 
polynomially with the number of spins in the lattice. However, the precision with which the Fourier 
components must be calculated is found to be an exponential function of the system size. 
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I. INTRODUCTION 



Eq. ([!]) may be written in the form 



Spin lattice models are useful for the study of magnetic 
ordering in real materials. The dynamics of these models 
are specified by a Hamiltonian TL involving spin opera- 
tors for each of the n lattice sites. Of particular interest 
is the behavior of thermodynamic functions - such as the 
magnetization, specific heat capacity, and magnetic sus- 
ceptibility - across phase transitions. These functions are 
encapsulated in the dependence of the Helmholtz free en- 
ergy per spin, F, on system parameters such as the tem- 
perature or applied magnetic field; partial derivatives of 
F yield the thermodynamic functions. Thus, calculation 
of F over a wide parameter space suffices for the deter- 
mination of the finite temperature behavior of the spin 
lattice model. 

Calculation of the free energy for a general spin lattice 
by conventional means is difficult. A naive approach is 
to enumerate the eigenenergies {E m } of H, since 
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where k B = /3 _1 is the thermal energy, and Z is the 
partition function. However, as the number of eigen- 
states grows exponentially with the number of spins in 
the lattice, the time required to perform the calculation 
is prohibitively large. A variety of quantum Monte Carlo 
methods exist to calculate the free energy, including ther- 
modynamic integration G, 0K histogram methods |3[ [|] 
and cumulant expansion |5[ Q techniques. However, the 
"sign problem" (see, for example, Ref. |?]) prevents appli- 
cation of these methods to arbitrary lattice Hamiltonians. 

An alternate approach is available if one can efficiently 
generate an estimate of the density of states p(E). As 



F = -n- l k B 0\n 
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an approximation for the density of states p(E) directly 
translates into an estimate F for the free energy per spin. 

Algorithms for quantum computers have been pro- 
posed to determine information about the spectra of Her- 
mitian operators [§, g 0, 0, [TJ, |l§]. We focus on al- 
gorithms jl2], [l3| that efficiently generate estimates of 
individual Fourier components ft of p{E); they will be 
reviewed in detail in Section H. N iterations of the algo- 
rithms yield N Fourier components, from which an esti- 
mate of the density of states can be calculated. 

An important issue that has not been addressed is the 
efficiency of these algorithms for calculating thermody- 
namic functions as a function of n. For the calculation to 
be deemed efficient, it must be shown that the computa- 
tion time - and, thus, the number of Fourier components 
- required to calculate an estimate F to bounded error 
scales polynomially with n. The bounded error criterion 
we adopt is 



Prob (\F -F\< jk B 9 
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where 7 and e are small constants. Thus, the absolute er- 
ror in the estimated free energy per spin must be smaller 
than a fraction of the thermal energy with probability 
arbitrarily close to one. 

We examine the primary sources of error involved in 
the calculation of F to determine the efficiency of the 
spectral algorithms. First, as only a finite number of 
Fourier components fe of the density of states are calcu- 
lated, the estimated density of states is broadened rela- 
tive to the actual function. This deterministic source of 
error (i.e., it is unchanged if the calculation is repeated) 
is reduced by increasing the number of components N, 
and thus the computation time. Second, there is an in- 
herent stochastic source of error reflected in deviations 
in the estimated fe from the actual values. This error 
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could arise from imprecise implementation of logic gates 
or noise in the measurement process. 

In this paper, we will show that if all of the fa are 
known exactly, the bound in Eq. (||) may be met by a 
number of Fourier components that scales polynomially 
with n. Thus, the error due to the broadening of the den- 
sity of states does not prevent efficient estimation of the 
free energy per spin. However, the free energy becomes 
increasingly sensitive to random errors in each of the fa 
as the number of spins is increased. We will show that the 
precision of the output of the quantum algorithm must 
improve exponentially with n in order to sustain the con- 
dition in Eq. (|). Thus, it is questionable as to whether 
spectral algorithms can be applied to the calculation of 
thermodynamic functions. 

The paper is organized as follows: Section II reviews 
the quantum algorithms used to generate the components 
and discusses assumptions and expected properties of 
the spin Hamiltonian. Section III describes the calcula- 
tion of F from the Fourier components, and discusses the 
influence of sampling and window functions. In Section 
IV, we analyze the deterministic error due to broadening 
and determine the number of samples required to meet 
Eq. (^). In Section V, we analyze the impact of random 
deviations in the components fa on the estimated free 
energy. 



II. QUANTUM ALGORITHMS 

In this section, we review quantum algorithms for the 
calculation of the Fourier transform of the density of 
states. We describe a simple algorithm applicable only 
to Hamiltonians that are diagonal in the computational 
basis, and then discuss a more general algorithm p4| ap- 
plicable to ensemble quantum computers. 

In regards to notation, we use the standard model 
for quantum computation, assuming our p qubits to 
be two-level systems with logical states |0) ■ and 
j €E {0, 1, . . . ,p— 1}, corresponding to eigenstates of the 
a^P Pauli spin operator with eigenvalues ±1. The com- 
putational basis states for the quantum computer are de- 
noted as |x) = \x\) x 1^2)2 •■■ \ x p) p i where {xj} are the 
binary digits of the integer x. It is assumed that the 
quantum computer is capable of implementing a univer- 
sal set of elementary single-qubit and two-qubit gates. 
The evolution time of these gates is an implementation- 
dependent constant, such that the overall computation 
time is reflected by the number of gates used in the al- 
gorithm. 

We will restrict our discussion to lattices of spin- 1/2 
particles, as it is straightforward to map the eigenstates 
of in the spin system to the logical |0) ■ and states 
of qubit j of the quantum computer. It should be noted 
that this restriction does not preclude the treatment of 
lattices of particles with spins larger than 1/2. Gener- 
alized Jordan- Wigner transformations exist jl5|, |l6| to 



represent the dynamics of such lattices by a collection of 
spin-1/2 particles via intermediate fermionization. 

Prior to the discussion of individual algorithms, we 
state three assumptions regarding the nature of the spin 
lattice Hamiltonian. First, we assume that the en- 
ergy bandwidth AE — the energy difference between the 
ground state and the highest excited state - is bounded 
by a polynomial function of n. This assumption is likely 
to be valid for models of interest. As an example, con- 
sider a lattice of particles interacting by an nearest- 
neighbor XXZ interaction: 



H = 



(id) 



(4) 



The expectation value of the summand in Eq. (4|) must 
lie between -(2|J X | + \ J Z \) and (2|J X | + | J z \). The num- 
ber of terms in the summation is equal to n/2 times the 
coordination number for the lattice, implying that the 
energy bandwidth is bounded by a function linear in n. 
In general, for any Hamiltonian involving only pairwisc 
interactions [p2[ between spins (of n-independent inter- 
action energy), it is evident that the energy bandwidth 
is 0{n 2 ). 

Second, we assume that the time-evolution operator 
U(t) = exp(— iH£) can be implemented as a sequence of 
elementary single-qubit and two-qubit gates, where the 
number of gates is a polynomial function of n. In cases 
where the Hamiltonian consists of commuting pair- wise 
interactions {e.g., the Ising model), this decomposition is 
trivial. Otherwise, one may use a Trotter-Suzuki expan- 
sion of non-commuting terms 0, [TJj to implement U(t) 
to arbitrarily small error. 

Finally, we assume that the energy scale is shifted such 
that the eigenenergies fall between E = and E = AE. 
This last assumption is made for mathematical conve- 
nience, and does not affect the results of our analysis. 

The following algorithms are based on the fact that the 
Fourier transform of the density of states p{E) is equal 
to the trace of the time evolution operator: 



/(*) 
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As |/(£)| < 2™, it is convenient to define a function 



9(t) = ^ /(*) 



-Tr 



(5) 
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such that \g(t)\ < 1. The algorithms described in this 
section calculate samples of g(t) at discrete times fa. 

Before discussing the general case, it is illuminating 
to examine a simple algorithm restricted to spin lattices 
for which TL is diagonal in the computational basis. As 
an example, one might consider a nearest-neighbor Ising 
model in a longitudinal magnetic field: 
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FIG. 1: Logic diagram of an elementary algorithm to estimate 
Iff* I- 



The gates shown in Fig. 1 for an n-qubit computer can 
be used to calculate the magnitude of g(tt) = ge- The 
quantum computer is initialized to the |0) state. The 
gate W corresponds to a sequence of Walsh-Hadamard 
gates Wj for each qubit j: 
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The gate U(te) corresponds to the time evolution opera- 
tor. 

As a final step, a projective measurement is performed 
in the computational basis. It is straightforward to show 
that the probability of observing all qubits in the logical 
|0) state is equal to \ge\ 2 - 
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By assumption, the computational basis states \m) are 
eigenstates of the Hamiltonian, and the time-evolution 
operator appends a phase proportional to the eigenvalue 
E m to each term. An unbiased estimator for \gi\ can be 
derived by performing R repetitions of the algorithm R, 
and counting the number of times all qubits are found in 
the logical |0) state. 

The magnitude of gi is insufficient to reconstruct the 
density of states. By adding an ancilla qubit a, as shown 
in Fig. 2, one may extract estimates of both the real 

and imaginary parts of ge. The X = exp (iTta^ /4j 

and Y = exp (iira^ /4j gates correspond to tt/2 rota- 
tions of the ancilla qubit. If the X gate is used, the 
probabilities of observing the |<fo) = |0) a |0) g . . . |0) e 

l^> = |i>jo> gi " 



or 



|0) states are 



1 + ige 


2 


1 - ige 


2 


, PX1 = 


2 



Pxo 



respectively. The Y gate leads to probabilities 





1 + 


2 


i - ge 


PYO = 


2 


) PY1 = 


2 



(10) 



(11) 



a ■ 

Q2- 
Q3- 

qn- 



W„ 



W 



X/Y- 



u{ti) 



W 



FIG. 2: Logic diagram of an algorithm to estimate the real 
and imaginary parts of gi for diagonal 7i. 



By executing R iterations of the algorithm with the X 
gate and R iterations with the Y gate, one can derive 
estimators p for the probabilities. Unbiased estimates of 
the real and imaginary parts of gi are 
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We use the tilde to distinguish estimates of the Fourier 
components obtained from the quantum algorithm from 
the exact values. 

The algorithm described above depends on our ability 
to construct an equally-weighted coherent superposition 
of eigenstates of H; hence, the restriction to Hamiltoni- 
ans that are diagonal in the computational basis. It can 
be generalized jl4j by considering an algorithm involv- 
ing an ensemble of quantum computers, such that the 
ancilla qubit is still initialized to the |0) state, but the 
remaining n qubits are in a completely random mixed 
state. The initial density matrix for the system is 



i T fiW + ^)MW...jW 1 
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where 1^ is the identity operator for qubit I. The op- 
erator 

7(91)7(92) . ../(<?») is equal to the resolution of the 
identity J2e where {|^}} is an orthogonal set of 

states in the subspace spanned by qubits q± through q n . 
One could use as the eigenstates of the Hamilto- 

nian ft. We do not need to explicitly solve for the eigen- 
states of 7i; the initial density matrix can be thought of 
as an incoherent mixture of eigenstates for any Hamilto- 
nian TL. 

If the coherent superposition created by the Walsh- 
Hadamard gates is replaced by such an incoherent mix- 
ture, then an algorithm nearly identical to the one shown 
above works for any choice of 7i, as shown in Fig. 3. The 
final measurement is the expected value of averaged 
over the ensemble. 

If the X gate is used for the ancilla qubit, the expected 
value of af^ is 
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FIG. 3: Logic diagram of an ensemble algorithm to determine 
ge for arbitrary 7i. 



whereas the Y gate leads to 



Re («?,). 



(16) 



Thus, the real and imaginary components of the esti- 
mator gi can be calculated from two iterations of the 
ensemble quantum algorithm. 

The ensemble algorithm is attractive for two reasons. 
First, it is applicable to any spin-1/2 lattice Hamilto- 
nian H, provided that the time-evolution operator can 
be decomposed into a sufficiently small number of ele- 
mentary gates. Second, initial state preparation lends it- 
self to ensemble quantum computation proposals involv- 
ing spin resonance, where the Zeeman splitting between 
qubit spin states is small compared to the thermal energy. 
In equilibrium, the initial density matrix of the system is 
well-approximated by the identity operator. The single 
pseudopure state qubit can be created from two thermal 
spins |f§. 



III. FREE ENERGY ESTIMATION 

In this section, we discuss how an estimate F of the 
free energy is generated from the Fourier components of 
the density of states, and we examine the effects of dis- 
cretization on the estimated density of states. 

A conceptual overview of the free energy calculation 
including post-processing is shown in Fig. 4. N samples 
of f(t) are estimated via the quantum algorithm p3| , and 
are weighted by a windowing function b&(t), described 
below. Fourier transformation yields an estimate p{E) 
for the density of states. The density of states may be 
integrated to compute the partition function, and, thus, 
the free energy. 

As iterations of the quantum algorithm yield discrete 
samples of f(t), the reconstructed estimate p(E) is dis- 
torted relative to the exact function p(E) . This distortion 
translates into error in F. It is convenient to view this er- 
ror in the context of windowing and sampling of the exact 
Fourier transform f(t) of the density of states. Trunca- 
tion of f(t) to a window of width To centered about t = 
(i.e., multiplication of f(t) by a windowing function b\(t) 
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FIG. 4: Block diagram of the calculation of F. 



that is constant for \t\ < To/ 2 and zero elsewhere) leads 
to a convolution of the density of states by a broadening 
function b x (E) = cosine (|§) = a x [sin (^)] / (^), 
where the energy resolution Ae is given by 



Ae = 



2tt 



(17) 



The window is scaled such that the broadening function 
is normalized to unit area; i.e., ot\ — 1/Ae. Increasing 
the window size To reduces the width of the broadening 
function, and thus the error in the estimate of F . 

The effect of sampling on the estimated density of 
states can be determined by multiplying f(t)b\(t) by an 
impulse train s(t) of spacing At: 



s(t) = At S{t~£At). 



(18) 



Sampling leads to periodic replication of the broadened 
p(E). The resultant density of states is given by the 
Fourier transform of /(i)&i(t)s(i): 
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To avoid aliasing in the estimated density of states, the 
Nyquist sampling condition requires that 



At < 



2tt 
AE' 



(20) 



The spacing between samples of fit) is determined solely 
by the estimate of the energy bandwidth, AE. We as- 
sume that sampling is performed at the Nyquist rate, and 
the equality holds in Eq. (p0|). 

As the number of samples is equal to the ratio of the 
windowing function width To to the sampling time, one 
could determine the minimum value of T required to 
satisfy Eq. (||) as a function of n. However, the rect- 
angular windowing function leads to poor results. The 
envelope of the associated broadening function bi(E) falls 
off weakly as 1/E; the oscillating side lobes are amplified 
at low energies by the Boltzmann factor in the calcu- 
lation of the free energy. The window width required 
to mitigate the resultant error scales poorly with n. In 
contrast to using wider rectangular windows, one may 
adopt more elaborate window shapes whose correspond- 
ing broadening functions exhibit envelopes that are more 
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sharply peaked. We consider the functions b&(t) formed 
by the successive convolution of rectangular windows, 
each of width Tq. O is referred to as the order of the 
windowing function. For = 2, the window is triangu- 
lar and of width 2Tq. With increasing order, the window 
approaches a Gaussian shape, and is of width OXo- The 
resulting broadening function is then 



b e (E) = ae 



1 at-J 



(21) 



which exhibits a 1/E e envelope. The value of ae is 
determined, as the area under be(E) is one. In practice, a 
given window shape is constructed by obtaining samples 
ft within the window width QTq centered at t = 0, and 
weighting each sample by b@ t t = &e(t^). 

As the envelope of the side lobes of be (E) falls off ex- 
ponentially with 0, windowing functions of large order 
significantly reduce the error in the calculated free en- 
ergy. However, the tradeoff is a wider window, leading 
to more Fourier components, and thus more iterations of 
the quantum algorithm: 



N = 



QT 

At ■ 



(22) 



Therefore, the question of how N scales with the num- 
ber of spins, n, translates into the determination of the 
minimum values of 6 and To required to satisfy Eq. (||). 

An estimate Z for the partition function can be cal- 
culated directly from the estimated Fourier components 
without intermediate calculation of the density of states. 
We denote quantities obtained from the quantum algo- 
rithm with a tilde, in contrast to their exact values. First, 
note that the Fourier transform of f(t)be(t)s(t) may be 
evaluated explicitly via Eq. ([l8]) to give an estimate of 
the broadened, periodically replicated density of states 
p'(E) in terms of the components ff. 
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where we have defined tg = £At, and the sum is per- 
formed over all tt within the window described by Red- 
integrating Eq. (|23|) over the energy bandwidth and us- 
ing Eqs. (ra) and (Ed), one finds 
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Note that an estimate F of the free energy may be ob- 
tained from the logarithm of Eq. ( |24| ) . 

In addition to describing how an estimate of the free 
energy per spin is calculated from the Fourier compo- 
nents, Eq. (|24|) will serve as a starting point to determine 
the probabilistic error in F due to imprecise values of gt . 



IV. ERROR ANALYSIS: BROADENING 

In this section, we determine an upper bound on the 
number of samples N of g(t) required to calculate the 
free energy to the tolerance prescribed by Eq. (||). At 
this point, we consider the individual samples of j]£ to 
be known exactly, and only consider the error in F due 
to the finite number of Fourier components - i.e., due 
to broadening of the density of states. With this restric- 
tion, we can show that N is a polynomial function of the 
number of spins n. 

As it is more convenient to work with the partition 
function than the free energy, we use a more stringent 
bound based upon the relative error in the calculated 



partition function Z . As 



F — F 



e -7» _ 1 < _ _ < e 7« 



it is sufficient to demand that 
Z-Z 



Prob r = 



<n>i 



(25) 



(26) 



where £ = 1 — exp(— jn). Note that £ < 1, and in the 
limit 771 <C 1, £ — > 771. 

By Eqs. @, © and (||), if the Nyquist sampling 
condition is satisfied, then 



N 



®AE 
Ae ■ 



(27) 



It has been asserted that AE is a polynomial function 
of n. In the remainder of this section, we examine the 
dependence of O and Ae on n such that Eq. ( p6| ) is 
satisfied. We require a pair of intermediate results: 



Lemma 1: If be{E) (as defined in Eq. (pl|)) is subject 
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to the normalization condition 1 = b@(E)dE, then 



or / 
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(28) 



where c w 2.0367. 



Lemma 2: 



/Ae 
b&(E)dE<- 
-Ae * 



e-3 



(29) 



where is an even integer. 

Lemma 1 places an upper bound on ae such that 
b@(E) is normalized. Lemma 2 states an upper bound on 
the area ofbe(E) that is outside of the interval [— Ae, Ae] 
(i.e., outside the main lobe of the broadening function) 
that decreases exponentially with 0. Both are proved in 
the Appendix. Note that the proof of Lemma 2 applies 
to the case where O is even. 

One can relate the relative error r in the calculated par- 
tition function to the parameters and Ae via Lemma 
2. As the exact density of states p(E) may be expressed 
as a sum of delta functions for each eigenenergy E m , Eqs. 
(|l9|)-(^0|) at the Nyquist condition yield 

rAE 

Z = p{E)e~l 3E dE 
Jo 

^ oo pAE 

= Yj / b e (E - E m + kAE)e~ ,3E dE 

m u — — 
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(30) 



The change of variables allows one to view Z m as an in- 
tegral of the broadening function, centered at E m , and 
weighted by periodically-replicated segments of an expo- 
nential function. Z is found by summing over all eigenen- 
ergies. 

The maximum relative error r in the partition func- 
tion is bounded by the largest contribution r m = 



n i 



max 



z m -z„ 



from any single eigenenergy, where Z„ 



Define 7™ = Z m /Z m . Then, 
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IEm(7m - 1)Z„ 



< max|7 m - 1| = 



Em Z m 

Z m ~ Z n 



Z n 



(31) 



This argument shows that one may consider a sim- 
plified system with just one eigenstate at an energy E m 
somewhere in the energy bandwidth. An upper bound on 
the error r for this simplified system for any E m suffices 



to bounded the error for an arbitrary energy spectrum 
over the same bandwidth. 

Lower and upper bounds on Z m (Z m ,m\n and Z mtiaa , x , 
respectively) are now derived to bound r m , since 



< max 



7 — 7 



(32) 



In the main lobe, the minimum value of the Boltzmann 
factor is e -/ 3 ( s ™+ Ae ). Outside of the main lobe, the min- 
imum value is e~@ AE . Thus, 
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(33) 



Similarly, as the maximum value of the Boltzmann factor 
is e -?{ E m-Ae) jjjgide ^ ne ma i n i be and one outside, 

Z m < (1 - A sidc )e- f3 ^- Ae ^ + ^0 = Z m , max . (34) 

Substituting Eqs. @ and (§|) into Eq. ©, we see 
that 



r m < max[l - (1 - A sidc )e 0Ae - A sido e" 



-j3(AE-E m ) 



(l-A sMc )^ Ae + A sidc e^-l]. 



(35) 



It is difficult to invert Eq. ( J35| ) explicitly to find optimal 
conditions on A S id c (©) and Ae that ensure that r m < £. 
However, one can show that the following conditions are 
sufficient: 



(36) 
(37) 



(38) 



/3Ae = ln(l+£/2), 
A sido < \e^ AE . 

As proof of their sufficiency, note that 

l-(l-A sidc )e-^~A sidc e-^ AE - E - 

and, 

(l-A side )e' 3Ae + A sidc e^ -1 

< 1 4. i p -P(AE-E m ) 

2 2 

< (39) 

Therefore, the conditions in Eqs. ( j36|) and (|37]) guarantee 
that r < r m < £, as desired. 

Using Lemma 2, one can manipulate Eqs. ( j36| ) and 
([37]) to show that A scales polynomially with n. 



A Ml+l/2) 
Q _ ln^ > ^ + ln M ) + ^ 



2 In 7r In 7r 



ln7r 



(40) 
(41) 
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where k = 5/2 + '"^//^ w 2.9443 . As In 6 < 6, a 
sufficient condition to satisfy Eq. ( pl| ) is 



0/2 = \nfiAE + fj, hx(l/0 + k'] , 



(42) 



where // = 1/(2 In 7T — 1) and k! = [iKhnr. 

In summary, the error bound on the partition function 
is satisfied if the energy resolution scales linearly with 
temperature, and if scales linearly with f3AE. 

As a final step, we substitute the conditions in Eqs. 
@ and @ into Eq. @, disregarding the weak loga- 
rithmic dependence of and Ae on n: 



QAE 
N = — — oc 
Ae 



(J3AE)(AE) 

VP 



I3 2 (AE) 2 ex poly(n), 



(43) 



by the assertion that the energy bandwidth is a polyno- 
mial function of the number of spins in our system. This 
result shows that in the absence of error in the calculated 
Fourier components of the density of states, the free en- 
ergy per spin can be determined efficiently to bounded 
error. 



V. ERROR ANALYSIS: FOURIER 
COMPONENTS 

We now consider random errors in the individual val- 
ues of t/£, which may arise from imprecise implementa- 
tion of logic gates, or noise in the measurement process. 
Treating Re(<^) and Im(^) as random variables, these 
fluctuations are modelled by their variances. We assume 
that the variances a 2 are independent |24[ of I. In this 
section, the dependence of the maximum allowable value 
of a 2 on n such that Eq. (|^) is maintained is derived. 

As the estimate for the partition function Z is a linear 
combination of the independent random variables Re(g^) 
and Im(^), the variance of Z can be calculated from Eq. 



4=4" 



/ At \ 2 ^ 



(1-e- 



9^1 



(44) 



than 5 x 10 3 for > 40. Second, it is assumed that 
PI At = f3AE/2n > 1. The energy bandwidth is thus 
much larger than the thermal energy. This condition 
assures that the sum can be well-approximated by the 
integral 



\W) 9 J 1 



dt 



t 2 /(3 2 At 



4V 



(3AE 



(46) 



Eq. ( |4^ ) indicates that the standard deviation of Z 
scales exponentially (271] with n; i.e., as 2". Note that the 
exact partition function Z will typically be a more slowly 
increasing function of n. If the energy eigenvalues are 
limited to the domain [0, AE), then 2™ is an upper bound 
for the value of the partition function (achieved at infinite 
temperature, or if all eigenstates are degenerate with zero 
energy). Consider two simple examples. For the case of 
n non-interacting spins in a magnetic field with Zeeman 
energy h, Z = (1 + e ~P h ) n < 2"; for a linear chain Ising 
model in zero magnetic field, described by Eq. (0), Z = 
(1 + e - 2 J y i for periodic boundary conditions. Thus, if 
the distribution function for Z is Gaussian, one expects 
that the standard deviation increases exponentially faster 
|p8| than the mean Z. 

The above result may be used to derive a condition on 
a 2 such that the error bound on the free energy per spin 

is fulfilled. By Eq. (0), the condition \F - F\ < jk B G in 
Eq. (0) is equivalent to 



Ze" 7 " < Z < 2e 7 ". 



(47) 



Assuming a Gaussian distribution for Z centered about 
Z, 



e = l- Prob(Ze" 711 < Z < Ze in ) 



(48) 




Z(e 7 ™ - 1) 
V2a? 



erfc 



Z(l 



This result can be simplified if we consider the limit <§; 
1, such that e ±7 " s» 1±7«; i.e., for small desired absolute 
error in the free energy relative to the number of spins: 



If we assume that Z is Gaussian-distributed J25J, then 
the probability of Z deviating from its exact value Z can 
be related to the variance. Thus, the sum in Eq. (Q) is 
evaluated by making two simplifications. First, we model 
the windowing function b@(t) as a Gaussian. Recall that 
6e(i) is constructed by the convolution of rectangular 
windows of width To. In the limit of large 0, &e(^) may 
be approximated |26| by 



be(t) 



-t 2 /2u 2 



(45) 



where v 2 = 0Tq/12. Although this approximation over- 
estimates &e(0 away from t = 0, the fractional er- 
ror in Eq. (H) incurred by the approximation is less 



c = erfc 



The argument of the erfc(-) function must be order unity 
or larger for e < 0.1, so 



., = Q (Z 2 yo\y{n) 



(50) 



By the above argument, the variance in the measured 
Fourier components must decrease exponentially with n. 

The condition on a 2 is likely to translate into an expo- 
nentially scaling computation time for the overall calcu- 
lation. For example, consider as a quantum computer an 
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ensemble of spin-1/2 nuclei, where readout is performed 
by measuring the voltage induced in a pickup coil by free 
induction. A source of error in the measured Fourier 
components is the Johnson-Nyquist voltage noise due to 
the resistance of the coil [EG]. The variance in the ob- 
served voltage - and thus in the estimates for Re(ge) and 
Im(gi) - is inversely proportional to the measurement 
time. Thus, Eq. (|50|) implies that an exponentially long 
measurement time is required to satisfy the condition in 
Eq. (|). 

VI. CONCLUSION 

We examined the applicability of spectral quantum al- 
gorithms for the calculation of the free energy of spin 
lattice models. Provided that the time-evolution opera- 
tor for the system is decomposable into an efficient num- 
ber of elementary gates, an ensemble quantum algorithm 
exists to generate estimates of the density of states by 
calculating individual Fourier components of p(E). We 
analyzed the efficiency of this algorithm in calculating the 
free energy per spin of the system to bounded absolute 
error. 



The error in the calculated free energy arises from the 
calculation of only a discrete number of Fourier compo- 
nents fi, as well as from deviations in the measured val- 
ues of fi due to statistical errors. The first source of 
error, attributable to broadening in the estimated den- 
sity of states, was shown to lead to bounded error with a 
number of Fourier components that is polynomial in n. 
Thus, if the components fe are known exactly, the spec- 
tral algorithm is an efficient means to calculate the free 
energy per spin. However, the effect of random devia- 
tions in the calculated values of fi grows with increasing 
n. As the size of the system increases, the maximum 
tolerable variance in measured Fourier components de- 
creases as Z 2 /4™ for large n and small absolute error. As 
an upper bound for the partition function is 2 n , the spec- 
tral algorithms are not an efficient method to determine 
F in the presence of statistical errors in fi. 
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APPENDIX A: PROOFS OF LEMMAS 1 AND 2 

Proof of Lemma 1: A lower bound is first derived for 

/OO /"00 
[smc(x)] e dx = e eln[sinc{x)] dx. (Af) 
-OO J —OO 

We exclude infinitesimal regions around x — mir (m € Z) 
from the integral to avoid divergence of the logarithm; as 
sinc(a;) approaches a finite value in these regions, the 
contribution of these regions to the integral can be made 
arbitrarily small. 

Using a series expansion for In [sinc(x)] pjj, 



In [sinc(x)] 



00 x 2k / 00 ]_ ^ 
fc=2 \ra=l / 



> 



_.2A- 



6 / ^— ' kir 2k 

k=1 



Thus, 



7> / e - e * 2 /6 e - 



" s fc= 



k = 2 i-^Ik 



dx. 



(A2) 



(A3) 



The integrand is positive over the entire domain of a;, and 
both exponential factors monotonically decrease with \x\. 
Thus, one may place a lower bound on I by reducing 
the limits of integration to any finite interval, such as 
\x\< v/67e. Thus, 



I > 



\r~- oc (6/e) fc 
" 2_, fc = 2 ZZZh 



/6/e 



/6/e 



(A4) 



The integral is yfhr/G) erf(l). The summation can be 
performed explicitly to yield 

j > e i+^ei„(i-6A 2 e)/ 6 ^ crf(1) 



= e 1- 



97T 2 



en- 2 /6 



erf(l) 1 



> e 1 - 



erf(l) 



(A5) 



where we make use of the fact that (1 — l/x) x is a mono- 
tonically increasing function for x > 1. 

This lower bound for / is used to establish an upper 
bound for a e . 



Ae I 



< Ae \ V 67r 



(A6) 



where c is defined as 
f 



f 

e 



1 -6/vr 2 



77 2 /6 



erf(l) 



2.0367. 



(A7) 



Proof of Lemma 2: For even, bo(E) is a non- negative 
function with unit area. If one treats bg>(E) as a proba- 
bility density function, one can use the Markov inequality 
to bound the area outside of the main lobe. 

Consider a random variable Y with support y > 0; i.e., 
Y only takes non- negative values. Markov's inequality 
bounds the probability of deviations from the mean: 



Pr (Y > 5) < 



E(Y) 



(A8) 



where E(Y) is the expectation value of Y. Define a sec- 
ond random variable X, such that Y = [X — E(X)] m , 
where m is an even integer. Then, 



Pr{[X - E(X)} m >5}< 
=>Pr{\X -E(X)\ >e}< 



E {[X - E{X)] m } 
5 

E{[X - E(X)] m } 



(A9) 



This bound is expressed in terms of the m th central mo- 
ment of X, if it exists. The result reduces to Chebyshev's 
inequality for m = 2. 

Note that if one treats b&{E) as a probability distri- 
bution function for a zero-mean random variable E, the 
above inequality provides a bound for the area outside 
the main lobe (i.e., e = Ae). The central moment is 
evaluated for m = — 2. 

/OO 
E e ~ 2 b e {E)dE 
-oo 



< a e — 



Ae 
Ae 
Ae 

7T 



e-i 



e-i 



e-i 



sin x 



-dx 



sin 2 x 



-dx 



7T 



If we define the area outside the main lobe as 



^4-sidc — 1 



b e {E)dE, 



Ae 



then. 



A 



side 



Pv{\X - E(X)\ > Ae} < 



ae Ae 

■ 



Combining Eq. [A.12] with Lemma 1, we find 



r e-3 



(AI0) 



(All) 



(A12) 



(A13) 



